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I discuss the accuracy requirements on numerical relativity calculations of inspiraling compact 
object binaries whose extracted gravitational waveforms are to be used as templates for matched 
filtering signal extraction and physical parameter estimation in modern interferometric gravitational 
wave detectors. Using a post-Newtonian point particle model for the pre-merger phase of the binary 
inspiral, I calculate the maximum allowable errors for the mass and relative velocity and positions 
of the binary during numerical simulations of the binary inspiral. These maximum allowable er- 
rors are compared to the errors of state-of-the-art numerical simulations of multiple-orbit binary 
neutron star calculations in full general relativity, and are found to be smaller by several orders 
of magnitude. A post-Newtonian model for the error of these numerical simulations suggests that 
adaptive mesh refinement coupled with second order accurate finite difference codes will not be able 
to robustly obtain the accuracy required for reliable gravitational wave extraction on Terabyte-scale 
computers. I conclude that higher order methods (higher order finite difference methods and/or 
spectral methods) combined with adaptive mesh refinement and/or multipatch technology will be 
needed for robustly accurate gravitational wave extraction from numerical relativity calculations of 
binary coalescence scenarios. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.40.Dg, 02.60.Cb 
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I. INTRODUCTION 

Studies in numerical relativity in the past decades have 
claimed at least partial motivation from the imminent 
direct detection of gravitational waves from both ground 
based (LIGO, GEO, TAMA, VIRGO) and space based 
(LISA) detectors. Theoretical calculations of gravita- 
tional waveforms from realistic astrophysical phenom- 
ena will be an essential ingredient in the extraction of 
characteristic information of gravitational wave sources 
(e.g., mass, spin, size, and composition of compact ob- 
jects) from the detected gravitational waves. In particu- 
lar, gravitational waves produced during the coalescence 
of binary compact objects (neutron stars and/or black 
holes) are strong candidates for direct detection, and thus 
it is precisely these systems that are of great interest 
to the numerical relativity community. Recent advances 
in numerical relativity, in particular with respect to the 
stability of binary black hole evolutions 0, Q> an d bi- 
nary neutron star evolutions [3- 0- El ■ make possible the 
calculation of gravitational waves from fully general and 
consistent numerical relativity simulations of binary co- 
alescences. 

However, numerical relativity simulations contain er- 
rors that arise from attempting to solve continuum dif- 
ferential equations (the Einstein field equations) on infi- 
nite domains (asymptotically flat spacetimes) with digi- 
tal computers of finite size and speed. Examples of these 
errors are truncation errors (e.g., due to the truncating 
of Taylor series expansions for finite difference methods 
or to the truncating of function expansions for spectral 
methods) and boundary errors (e.g., errors induced by 
the introduction of a computational domain in causal 
contact with the binary and/or the gravitational waves 



being emitted). The magnitude of these errors deter- 
mines the accuracy of the numerical relativity simulation 
(here, I do not include modeling errors in the determina- 
tion of the accuracy of a numerical relativity code, e.g., 
inaccurate equation of state for neutron star matter or 
astrophysically incorrect initial data for binary simula- 
tions). In this paper, I demonstrate, for the first time, 
a calculation for determining the accuracy required of 
inspiraling binary numerical relativity simulations in or- 
der that the characteristics of the extracted gravitational 
waveform represent the physics of the binary system to 
the experimental error level of the gravitational wave de- 
tector. Using the gravitational waveform accuracy crite- 
rion in ,7] , I calculate the sensitivity of the gravitational 
waveform to various parameters of the dynamical binary 
system (e.g., binary separation, angular velocity, mass) 
assuming a specific target sensitivity for the gravitational 
wave detector; forcing the errors in the same dynamical 
parameters within numerical relativity coalescence simu- 
lations to be smaller than these sensitivities will be one 
way of guaranteeing an extracted theoretical waveform 
accurate to the sensitivity level of the gravitational wave 
detector. Truncation and boundary errors in the orbital 
separation of the multiple-orbit binary neutron star sim- 
ulations in are calculated and are shown to be several 
orders of magnitude larger than the margin allowed for 
by the gravitational wave sensitivity calculation. Using 
a post-Newtonian model of the truncation and bound- 
ary errors in the binary neutron star numerical relativ- 
ity simulations, I estimate the computational resources 
required for accurate gravitational waveform generation 
from such simulations, and conclude that both higher 
order methods and mesh refinement or multipatch tech- 
nology will be required for robust, reliable, and accu- 
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rate gravitational waveform extraction from numerical 
relativity simulations of coalescing binary inspirals and 
mergers. 

The outline of the rest of the paper is as follows. In 
Section ^ the sensitivity of gravitational waveforms to 
various physical characteristics of the coalescing binary 
is calculated for binary black hole and binary neutron 
star systems. In Section IHII I compare these sensitivi- 
ties to errors in recent state-of-the-art binary simulations 
in numerical relativity, and find the errors to be orders 
of magnitude larger than the gravitational wave sensi- 
tivities. I conclude by discussing methods that may be 
helpful in reducing the error of numerical relativity calcu- 
lations of coalescing binaries to levels that would permit 
their extracted gravitational waveforms to be used with 
confidence as templates in modern interferometric grav- 
itational wave detectors. 



II. REQUIRED ACCURACY FOR NUMERICAL 
RELATIVITY SIMULATIONS OF BINARY 
INSPIRALS: SENSITIVITY OF 
GRAVITATIONAL WAVEFORMS TO PHYSICAL 
CHARACTERISTICS OF THE BINARY 

In order to calculate the sensitivity of gravitational 
waveforms to various physical characteristics of the coa- 
lescing binary, I calculate binary coalescence solutions to 
the post-Newtonian equations of motion for non-spinning 
point particles. This formalism has the advantage that 
accurate gravitational waveforms can be calculated in as- 
trophysically relevant binary coalescence scenarios. The 
major disadvantage is that at small binary separations, 
the point particle approximation breaks down due to fi- 
nite size effects; the gravitational waveforms obtained 
by solving the post-Newtonian equations of motion must 
therefore be truncated at the point where these effects be- 
come important. As a result, the sensitivities calculated 
here will only bound the true sensitivity, i.e. the sensitiv- 
ity of the entire gravitational wavetrain, through plunge, 
merger, and ringdown of the final merged object. In other 
words, due to the fact that the gravitational waveform is 
being truncated when finite size effects become impor- 
tant, the sensitivity of the entire physical gravitational 
waveform to variations in the physical characteristics of 
the system is being underestimated. Thus, the errors in 
numerical relativity simulations must be at least as small 
as the sensitivities calculated here. 



A. Post-Newtonian equations of motion 

The general relativistic equations of motion for non- 
spinning point particles in harmonic coordinates with po- 
sitions xi and 2?2, and masses mi and ma, can be written 
in a post-Newtonian expansion as 

d 2 x m Tn 

-772 = 2 " + ~2 [H-A-tPN + MPN + A 3PN H ) + 



rv(B 1PN + B 2PN + B 3PN H )] + 

8 777 TiT 

-pV^r— [rn(A 2 . 5 pN + A 3 . 5PN + A 4 . 5PN H )- 

5 r z r 

v(B 2 .5PN + B 3 . 5 pjy + B^.spjy +•••)], (1) 

where x = x 2 — x\ is the relative separation of the parti- 
cles, v — v 2 — v\ is the relative velocity between the par- 
ticles, r — \x\, h — x/r, m = mi + m2, r\ = mim 2 /m 2 , 
and r — dr/dt. The post-Newtonian expansion in Eq.^is 
carried out in powers of e ~ ra/r ~ v 2 (I set G = c = 1), 
where each power of e represents one post-Newtonian 
(PN) order in the series. The 1PN and 2PN terms are 
standardfe.g., see HUES)- The 2 ' 5PN H 113 03 and 
3.5PN [l(| have also been completely determined. The 
3.0PN terms have just recently been calculated 0] up to 
one gauge-dependent constant. Employing an energy and 
angular momentum balance technique, the 4.5PN terms 
have been fixed modulo 12 free "gauge" parameters 01 
(appendix B of demonstrates that these free param- 
eters have a negligible effect on inspiral dynamics). Once 
the initial relative positions and velocities of the parti- 
cles are given, Eq. ^ is then solved numerically for the 
time evolution of the binary system. Specifically, if the 
initial separation r and its time derivative r along with 
the initial relative angular position (f> and its time deriva- 
tive (j> are specified at time t — 0, then the equations of 
motion Eq. ^ specifies r(t) and 4>{t) for t > (I assume 
the binary orbits within the z = plane). 

I use the post-Newtonian formalism presented in [TH 
0,^3,01 t° calculate the polarization state h(t) = h+(i) 
of the gravitational radiation as a function of the motion 
of the binary. For definiteness, I fix both observation an- 
gles ($, 9) to be for the remainder of the paper (i.e., the 
binary, which is orbiting in the z — plane, is observed 
along the +z-axis). 

B. Waveform sensitivity to physical characteristics 
of the binary 

An inner product on the space of waveforms h(t) is 
defined as 

where hi(f) and h 2 (f ) are the Fourier transforms of the 
two waveforms h\(t) and h 2 (t), and Sh(f) is the one- 
sided power spectral density of the strain noise of the 
detector. For the calculations in the remainder of this 
paper, the model of the one-sided power spectral density 
of the strain noise for the advanced LIGO detector found 
in 19] is used for Sh(f), where the mass scale is set by 
assuming an equal mass binary with mi = m2 = 1.4 M@. 
The latest detection rates for binary neutron star coales- 
cences for the advanced LIGO sensitivity is between 40 
and 650 events per year [2(| • Assume that the gravita- 
tional waveform h(t) contains an error Sh(t). Arguments 



in Q determine the criteria that the error Sh(t) be small 
enough so that the quantity 

a _ 1 (Sh\Sh) 

A " 2 (h\h) [6) 

satisfies 

A < 0.01. (4) 

This accuracy criteria is based on matched filtering ac- 
curacy arguments, physical parameter estimation argu- 
ments, and arguments based on the total information 
content of the gravitational wave signal (see for de- 
tails). The specific criteria Eq.0]is based on the assump- 
tion that the signal to noise ratio is of order 10. Much 
higher signal to noise ratios are expected for LISA, and 
therefore an even more stringent requirement on Sh(t) 
could be required in that case (e.g., for gravitational wave 
templates used in physical parameter estimation, the re- 
quired A scales like the inverse square of the signal to 
noise ratio). 

Using the post-Newtonian equations of motion de- 
scribed in Sec. Ill Al as a model for coalescing binary 
compact objects, along with the gravitational wave ac- 
curacy requirements Eqs. |3| and I now analyze the 
sensitivity of the gravitational waveform to various phys- 
ical characteristics of the binary dynamics. Note that 
the gravitational waveform h(t) of our post-Newtonian 
model of a binary inspiral depends solely on the 8 pa- 
rameters A = {ro,<po,ro,(po,mi,m2,tf,R}. Here, ro and 
4>o specifies the relative binary position at the initial time 
t = to, ro and (f>o specifies the initial relative binary veloc- 
ity, mi and mi specifies the mass of each compact object, 
while R specifies the observation distance from the center 
of mass of the binary. The parameter tf is the time at 
which the waveform is truncated due to the breakdown of 
the post-Newtonian point particle approximation. This 
breakdown occurs when the internal structure of the in- 
dividual compact objects has a significant effect on the 
orbital dynamics of the binary. Throughout the rest of 
this paper, I take tf to be the time at which the binary 
separation r(tf) = 4m for equal mass binary black holes 
and r(tf) = 8 m for equal mass binary neutron stars. 
I orient the system such that the initial relative angle 
0o = 0. For a specific initial binary separation r , I set 
the initial relative velocity parameters ro and 4>o to be 
those specified by the unique quasicircular solution of the 
binary (the quasicircular solution to the PN equations of 
motion is obtained by starting from circular orbit initial 
data in the limit as the initial separation r — ► oo, see 
Sec. 2 of [3). 

For a specific set of parameters A that specify the wave- 
form h(t), define the quantity Ao.oiA to be the ranges in 
the parameters such that the change in the gravitational 
waveform Sh(t) induced by separately changing each in- 
dividual component A^ satisfies A < 0.01. The quan- 
tity Aq.oiA can be interpreted as the "allowable" error in 
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FIG. 1: An example of the definition of the Ao.oi operator. 
A (Eq. is plotted as a function of the variation of the ini- 
tial orbital separation Sr for an equal-mass black hole binary, 
starting 3 orbits before merger (m is the total mass of the sys- 
tem). In order that variations in the resulting gravitational 
waveform satisfy A < 0.01, any variation in the initial orbital 
separation r must satisfy Sr < 0.0165 m. I therefore define 
Ao.oir = 0.0165 m. 



each parameter during the course of a numerical relativ- 
ity time evolution simulation, since changes within this 
range do not appreciably (to the tolerance set by Eq. |3J) 
affect the gravitational waveform. 

As a concrete example, I examine the sensitivity of the 
gravitational waveform on changes in the binary separa- 
tion r. I solve the 4.5PN equations of motion for the qua- 
sicircular solution for an equal mass (mi = mz = m/2) 
binary and find that exactly 3 orbits before merger (for 
definiteness, assume binary black holes and define the 
merger to be at binary separation r — 4 m) , the binary 
separation is r = 7.3195 m, the relative radial velocity 
is r = —0.005573, and the relative angular velocity is 
4> = 0.04263/m. Using these values of parameters as ini- 
tial data, the subsequent gravitational waveform h(t) for 
the last 3 orbits before merger is computed. I define a sec- 
ond gravitational waveform h'(t) to be the one obtained 
by changing the initial binary separation ro — > ro + Sr, 
keeping all other initial parameters fixed. The change in 
the waveform Sh(t) induced by changing the initial binary 
separation an amount Sr is therefore Sh(t) = h'(t) — h(t). 
The quantity A (Eq. [21) can now be calculated; in Fig.^ 
the quantity A is plotted as a function of the change in 
initial binary separation Sr. The quantity Ao.oiT is de- 
fined as the range of Sr such that A < 0.01; in this case, 
Ao.oiT = 0.0165 m (see Fig.^). The intuitive interpreta- 
tion of this calculations is as follows: Sr represents the 
error in the binary separation of a numerical relativity 
simulation of equal mass black holes at a time when 3 
orbits remain until the merger. The quantity Ao.oiT rep- 
resents the "allowable" error in binary separation r as 
set by the tolerance level of Eq. 0] Thus, a numerical 
error \Sr\ in the binary separation of a numerical rela- 
tivity simulation at a time when 3 orbits remain until 
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FIG. 2: Gravitational wave sensitivities to various dynami- 
cal quantities for orbiting binary black holes (top panel) and 
binary neutron stars (lower panel) are plotted as a function 
of the number of orbits remaining until merger. Shown are 
gravitational wave sensitivities of the mass of each compact 
object (Ao.oiwii = Ao.oi^ia, solid line), the orbital separa- 
tion (Ao.oiT, dotted line), the relative orbital angular velocity 
(Ao.oi0, dashed line), and the relative orbital radial velocity 
(Ao.oiX, alternating dot-dashed line). 



III. CASE STUDY: BINARY NEUTRON STAR 
EVOLUTION 

A. comparing errors in numerical relativity 
simulations of coalescing binary neutron stars to 
gravitational waveform sensitivities 

Several recent studies 0, have analyzed various as- 
pects of orbiting binary neutron stars by performing 3+1 
general relativistic hydrodynamic simulations. It is in- 
structive to analyze the accuracy of these simulations in 
comparison to the sensitivity required for the accurate 
generation of gravitational waves as represented in Fig. [2] 

In Q, quasiequilibrium initial data sets correspond- 
ing to two equal mass, corotating binary neutron stars 
were numerically evolved using a general 3+1 numerical 
relativity code that simultaneously solved the Einstein 
field equations coupled to the general relativistic hydro- 
dynamics equations. In order to assess the quality of 
the numerical solution, the exact same initial data set 
was numerically evolved 5 separate times, using a variety 
of discretization parameters and outer boundary place- 
ments (see Table II in Q). There are two types of nu- 
merical errors associated with the simulations presented 
in p|. The first type of numerical error is due to the 
finite difference approximation, where derivatives in the 
partial differential equations have been replaced by trun- 
cated Taylor-series approximations. Assuming that the 
finite difference equations used in £| are consistent and 
stable, Lax convergence theorems (see [2lJ) state that the 
relationship between any quantity Q constructed from 
a true solution of the differential equation is related to 
that same quantity observed in the numerical solution 

Q numerical by 



= 



numerica 



i + d{Ax) + C 2 (Ax) 2 + 



(5) 



where Cj are constants (which are different for distinct 
quantities 0) and Ax is the discretization parameter used 
in the construction of the finite difference equations. As- 
suming that Ax is made small enough so that higher 
order terms can be ignored, define the truncation error 
{AQ) trunc of the calculations in to be 



(A0) t 



d(Ax) + C 2 (Ax) 



(6) 



merger that is greater than A0.01T would correspond to 
an unacceptable level of error in that simulation. 



In Fig. |21 the parameter tolerance An. 01 K for physical 
parameters mi , to 2 , rn , 4>q and f is plotted as a function 
of the number of orbits until merger for equal mass bi- 
nary black holes and binary neutron stars. Note that the 
gravitational waveform is roughly an order of magnitude 
more sensitive to small changes in the angular velocity 
as compared to small changes in the radial velocity. 



While the code being analyzed is formally second order 
convergent in both space and time, the use of high reso- 
lution shock capturing (HRSC) methods renders the hy- 
drodynamics convergence rate to be first order in space 
in regions where the dynamical variables obtain a local 
extrema (see, e.g., 22]). Thus, the form of the trunca- 
tion error, Eq.|Sl necessarily contains a term proportional 
to the first power of the spatial discretization parameter 
Ax. 

The second type of numerical error that exists in the 
simulations presented in Q| is due to the location of the 
boundary of the computational domain. Ideally, the com- 
putational domain would be placed many gravitational 
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FIG. 3: The absolute value of the truncation error (Ar), 

V Jtrunc 

and boundary error (Ar) botln<i of the orbital separation is 
plotted as a function of the number of orbits for the binary 
neutron star simulations in Q. For comparison, the gravi- 
tational wave sensitivity to variations in the orbital separa- 
tion, A e r is also plotted. Assuming the use of gravitational 
waveforms produced from numerical simulations as matched 
filtering search templates in a gravitational wave detector, 
sensitivities of Ao.oi, Ao.i, and A0.5 correspond to a detector 
event loss rate of 3%, 27%, and 87%, respectively. 



wavelengths away from center of mass of the orbiting 
binary, thus minimizing the effect of the boundary on 
the dynamics of the binary. However, computational re- 
source limitations coupled with the fact that the code 
used in |4| is only second order accurate and does not 
employ adaptive mesh refinement, prevented the location 
of the outer boundary of the computational domain to be 
placed no farther than 1/4 of a gravitational wavelength 
from the center of mass of the system. I model the error 
associated with the close proximity of the boundary of 
the computational domain analogously with that of the 
truncation error: assume that the error in any particu- 
lar quantity Q induced by the boundary goes to as the 
distance between the center of mass of the system and 
the location of the outer boundary (denote this distance 
as Tb) goes to infinity. I expand this boundary-induced 
error, (AQ) bound , as a power series about — 00 and as- 
sume that the values of rt used in the calculations in Q 
are large enough so that the power series can be trun- 
cated as 



(^Q)bound - n + r&2 



(7) 



The total numerical error, therefore, is represented as the 
relationship between the exact quantity Qexact specified 
from the solution to the differential equations and quan- 
tity Q numerical produced by numerical simulations: 



Qexact — Qnumerical + (AQ) It 



(AQ) 



boujid' 



(8) 



I now compare the errors in the binary separation r 
contained in the simulations of Q with the sensitivity of 
the gravitational waveform to changes in the binary sep- 
aration, Ao.niT. I use Eq. [S] to calculate the numerical 



errors in the calculation of the binary separation r(t) for 
the simulations presented in [j. At each time t, I as- 
sume that the numerically computed binary separation 
^numerical takes the form of Eq. 



1 'exact — ^numerical ~t~ (A^)^ runc ~\~ (Ar)i. 



(9) 



Using the results from the five simulations displayed in 
Fig. 17 in |j| (which are reproduced in Fig. 0] of this 
paper), the five unknowns in Eq. (r exact , C\, C2, Bi, 
and B2) are specified at each time t. The absolute val- 
ues of the truncation error (Ar) trimc and the boundary 
error (Ar) 6 d of the binary separation are plotted in 
Fig. |3| For comparison, the gravitational wave sensi- 
tivity to binary separation, A .oi?* (see Fig. 0), is also 
plotted. Important to notice is the magnitude of the er- 
rors (both truncation errors and boundary errors) in the 
binary separation of the numerical relativity calculations 
in as compared to the sensitivity of the gravitational 
waveform to the binary separation: the simulation errors 
in the binary separation are several orders of magnitude 
larger than the gravitational waveform sensitivity to vari- 
ations in the binary separation! Assuming that they are 
to be used as signal detection and parameter estimation 
templates, the gravitational waveforms extracted from 
these numerical simulations will contain errors that dom- 
inate the experimental errors in modern interfcrometric 
gravitational wave detectors. Steps must be taken to 
improve the accuracy of these simulations before gravi- 
tational waves extracted from them can be considered for 
use as templates in gravitational wave detectors. 



B. estimating the accuracy and computational 
resources required to extract sufficiently accurate 
waveforms from binary coalescence simulations 

In order to obtain information on the computational 
resources required to reduce the boundary and trunca- 
tion errors of simulations such as those in Q down to 
acceptable levels (i.e., such that the errors in the gravi- 
tational waveforms Sh induced by these numerical errors 
satisfy A < 0.01), these errors and their effect on gravita- 
tional waveforms are modeled using the post-Newtonian 
equations of motion for spinless point particles, Eq. ^ 
The same mass and initial orbital separation is used in 
the post-Newtonian model as was used in the numeri- 
cal relativity simulations. The initial angular velocity of 
the binary is determined by the circular orbit assumption 
(see [3), 

which is consistent with the initial data used 
for the numerical relativity simulations in . In order to 
model the effects of truncation and boundary errors with 
the post-Newtonian model, it is modified in the following 
way. First, the initial angular velocity of the binary is set 
to be a function of discretization Ax and outer boundary 
placement r as 



,Ax, 



,Ax, 



(}>o(Ax,rb) = 4>circular+ <3\ ( ) + &l{ ) 

m m 



where ^circular is the angular velocity determined by the 
circular orbit assumption. Second, the post-Newtonian 
equations of motion are modified so that the evolution of 
the angular momentum is given by 



dL\ _ f dL \ +(J (- Aa; ) + cr ( Aa: ) 
dt J \dt J PN m m 

,m s ,m 2 
+a 7 (-) + a s (-) 11 

n n 



In the ideal limit of infinite resolution (Ax — ► 0) and 
infinite distance from the center of mass of the system 
to the computational boundaries (r& — > oo), the solu- 
tions to these modified post-Newtonian equations reduce 
to the standard post-Newtonian inspiral solutions start- 
ing with circular orbit initial data. The constants <j\ 
through a% are to be chosen so as to best reproduce the 
effects of the truncation errors and boundary errors from 
the multiple-orbit simulations of binary neutron stars 
in p|. To this end, I define an "error function" x( a i)> 
which measures the difference in the orbital separation 
profile in time between the 5 numerical relativity simu- 
lations NS-A through NS-E from Q (which I denote as 
r^t), j = 1,2,3,4,5) and the solution to the modified 
post-Newtonian equations Eqs. EJ ^3 and ^2 (which I 
denote as r pn (i, Ax, ?*(,, (Xi)), as 



1 



Jo f dt(r vn (t, (Ax) s , {r b ) p <Ti) - i4(f) 



(12) 

where (Ax) ■ and {rb)j are the discretization and bound- 
ary placement parameters used to produce the numerical 
relativity result rj( r (i) for each of the five multiple-orbit 
numerical relativity simulations of binary neutron stars 
performed in [4j . The time integrations in Eq. 1121 use 
tf = 610m, which corresponds to roughly 2.5 orbits of the 
binary system. A global minimum of x(<r™ m ) = 0.093m 
is found numerically by varying the ai parameters. A 
comparison among the five numerical relativity simu- 
lations NS-A through NS-E from and the solutions 
of the modified post-Newtonian equations of motion 
Eqs. n EH and II II using the cr™ m parameters that mini- 
mize x( a i) is shown in Fig.^J 

Using the solutions to the modified post-Newtonian 
equations of motion corresponding to parameters <j™ n 
that minimize x( f7 i) as a model for the effects of the 
truncation and boundary errors in the full numerical 
relativity simulations of orbiting binary neutron stars 
in [4(, I am now able to gauge the effect these errors 
have on the resulting gravitational waveform. Specifi- 
cally, the goal is to find bounds on the discretization 
parameter Ax and outer boundary location parameter 
rf, such that the error in the produced gravitational 
waveform satisfies, e.g., Eq. 0] In order to calculate 
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FIG. 4: Top panel: plotted is the orbital separation as a 
function of time from the numerical relativity simulations of 
orbiting binary neutron stars in |4|. These simulations were 
performed with the same initial data, but with different dis- 
cretization parameters Aa; and computational domain bound- 
ary placement parameters tv Shown are simulations NS-A 
(solid line), NS-B (dotted line), NS-C (short dashed line), 
NS-D (long dashed line), and NS-E (alternating dot-dashed 
line), evolved to time t = 610 m, which corresponds to 2.5 
orbits. Bottom panel: the modified post-Newtonian model 



r pn (t, (Ax).,(r b ). 



is plotted for the five discretization 



and boundary parameters {(Aa;) ., (r&) .}, (j = 1,2,3,4,5), 
corresponding to the parameters used in numerical relativity 
simulations NS-A through NS-E, respectively. The modified 
post-Newtonian model robustly encapsulates the effects of the 
truncation and boundary errors within the full numerical rel- 
ativity simulations. 



A = (l/2){Sh\6h) / (h\h), I take the target gravitational 
waveform hit) to be that determined by the solution to 
the modified post-Newtonian equations of motion in the 
limit as Ax — » and r\, — > 00, which is just the waveform 
obtained from the ordinary post-Newtonian equations of 
motion assuming initial data corresponding to a circular 
orbit. This waveform we denote as ho(t). The "error" in 
the waveform 5h(t) induced by the truncation error and 
boundary error in the numerical simulations can then be 
calculated as 



Sh(t, Ax, r b ) = h(t, Ax, r b ) - h (t) 



(13) 



where h(t, Ax, r b ) is the waveform obtained by the modi- 
fied post-Newtonian equations of motion using discretiza- 
tion parameter Ax and boundary placement parameter 
Th (and, of course, using the Oi parameters that minimizes 
x(cri)). Fig. 0is a plot of the target waveform ho(t) and 
the waveform h(t, Ax NS -A; t^ns-a); which corresponds to 
the best numerical relativity simulation NS-A (the solid 
line in Fig. B}. 

Shown in Fig. is a contour plot of A (Eq. 01 as a 
function of Ax and r\,. Contours for A = 0.1, 0.01, and 
0.001 are shown. The "peninsula" -like shape of the con- 
tours in Fig. El are due to the slightly offsetting effect of 
the truncation and boundary errors in the numerical sim- 
ulations of larger discretization parameters Ax tend 



7 




t/m 



FIG. 5: Gravitational waveforms (multiplied by the distance 
to the source, R). The solid line waveform corresponds to 
ho(t), which is the "zero error" waveform corresponding to 
the modified post-Newtonian solution with discretization pa- 
rameter Ax = and boundary placement parameter rj = oo. 
The dashed line waveform is h(t, Ax N s-k, t"(>ns-a)> which cor- 
responds to the numerical relativity simulation NS-A (see 
Fig. 0J. The difference between these waveforms (Eq. I13P 
corresponds to A = 0.054, which is 5 times larger than our 
target accuracy (see Eqs. I^Jand^J. 



to artificially increase the rate at which angular momen- 
tum is lost from the binary, while closer outer bound- 
ary placements (smaller r&) tend to have the opposite ef- 
fect. For reference, the computational memory resources 
for the unigrid numerical relativity code used in Q is 
shown in Fig. El indicating Gigabyte (1024 3 bytes), Ter- 
abyte (1024 4 bytes), and Petabyte (1024 5 bytes) require- 
ments. Adaptive mesh refinement (AMR) allows the effi- 
cient minimization of errors induced by the boundary by 
permitting the placement of the boundary of the compu- 
tational domain farther from the coalescing binary for a 
fixed amount of computational resources (other methods 
could also reduce boundary errors, such as employing a 
Cauchy-characteristic matching code [23| or using a null- 
approaching slicing far from the center of mass of the 
binary The computational memory resources for 

an AMR version of the code is also shown in Fig. (I 
have assumed that the finest resolution grid is the size of 
the compact objects, that the grid at each level has the 
same computational volume as every other grid, and that 
the grids are uniformly nested with a refinement ratio of 
2). 

An optimistic reading of Fig. [5] implies that numeri- 
cal relativity simulations using a 10 Terabyte computer 
would be able to attain the target accuracy of A = 0.01. 
However, the reliability and robustness of such a calcu- 
lation would be highly questionable, due to the fact that 
the discretization and boundary placement would have 
to be fine-tuned to reach the tip of the A = 0.01 contour 
peninsula in Fig. [S] In order to obtain a robustly accu- 
rate simulation from which gravitational waveforms could 
be extracted with confidence, it will likely be necessary 
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FIG. 6: Contour plots of the gravitational wave accuracy pa- 
rameter A (Eq. [3J as a function of discretization parameter 
Ax and distance from the center of mass of the binary to the 
computational domain boundary r& for the numerical relativ- 
ity simulations of binary neutron stars in y| . Configurations 
for various sized computers with a unigrid code and an AMR 
code (nested boxes, with grid refinement ratio of 2) is shown 
for reference. 



to, at a minimum, use a target resolution Aa; ta rget and 
boundary placement rb target such that the gravitational 
wave accuracy parameter A (Eq. |3J) satisfies A < 0.01 
for all Ace < Aa; target and r& > rf, target . From Fig. H3 we 
see that this minimum target configuration is at roughly 
Aa; ta rget ~ 0.002 m and rf, target ~ 2000 m. However, 
this minimum target configuration would not be possi- 
ble with a unigrid code, and would just barely be possi- 
ble with an AMR code on a Petabyte machine, although 
the execution time of such a simulation would render 
it highly impractical. In order to reduce the computa- 
tional resources required to perform sufficiently accurate 
inspiral calculations in numerical relativity, higher order 
methods will need to be employed in future calculations. 
Possible higher order extensions to the code in £| in- 
clude the use of spectral methods (where the truncation 
error drops off exponentially with the number of collo- 
cation points) or the use of higher order finite difference 
methods. Fig. [7| reproduces the results of Fig. [5] assum- 
ing that the truncation error of the simulation falls off 
as (Ax) 8 , which is consistent with using an eighth-order 
finite difference method. The increased accuracy of the 
eighth order method allows for a larger discretization pa- 
rameter Ax; the new minimum target discretization is 
Az targ et ~ 0.33 771. As seen in Fig. |7| a combination 
of both higher-order finite difference methods and AMR 
will yield a robustly accurate simulation on Gigabyte 
scale computers. Not only will the memory fingerprint 
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FIG. 7: Contour plots of A (Eq.0 as a function of discretiza- 
tion parameter Aa; and distance from the center of mass of the 
binary to the computational domain boundary r&, assuming 
an eighth-order finite difference method is used for the numer- 
ical relativity simulations of binary neutron stars. Configu- 
rations for various sized computers with a unigrid code and 
an AMR code is shown for reference. Note that higher order 
finite difference methods, couple with AMR, will be able to 
robustly obtain the required gravitational waveform accuracy 
A < 0.01 on Gigabyte computers. 



of a sufficiently accurate simulation be relatively small 
for an AMR code employing higher-order methods, but 
just as importantly, the execution time of a simulation 
using such a code will be considerably smaller than, e.g., 
the 200,000 CPU hours required for the NS-A simulation 
from [4J displayed in Fig. 0] 



C. generality of results 

Various aspects of the implementation of numerical rel- 
ativity simulations of coalescing binary compact objects, 
including details regarding the spacetime and hydrody- 
namics solvers, boundary conditions, initial data, gauge 
conditions, and total evolution times, could have a large 
impact on the details of the accuracy studies presented 
in sections IlII Al and lHI Bl For instance, the implementa- 
tion of constraint-preserving boundary conditions [25, 26] 
in the simulations presented in |4| could reduce by a 
significant amount the errors in the simulation induced 
by the outer boundaries. Of course, any approximation 
method employed in the generation of gravitational wave 
templates used for signal searches and parameter estima- 
tions in gravitational wave detectors must be validated; 
one must directly confirm that the approximation is good 
relative to the signal to noise ratio of the detector. As 
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FIG. 8: Top panel: the separation of binary neutron stars 
(normalized by the total baryonic rest mass of the system 
Mn) simulated by the general relativistic hydrodynamics code 
in |g| is plotted as a function of time for two orbital periods. 
For comparison, the solution to the post-Newtonian point par- 
ticle equations of motion (accurate to order (v/c) 9 ) for the 
same masses and initial conditions is shown. Bottom panel: 
the binary separation from the general relativistic simulations 
NS-A and NS-B in Q (also shown in Fig. fusing a different 
normalization) is plotted for two orbital periods. NS-B has a 
similar outer boundary placement as that used in the simula- 
tions of [fj, but NS-A has an outer boundary that is located 
twice as far away from the binary as that of simulation NS-B. 
The similarity of the top and bottom panels suggests that the 
unphysical increase in the binary separation of the "Stable" 
simulation in 6] is due to boundary errors, and that the mag- 
nitude of the errors are roughly the same for the simulations 
in H and M. 



such, the calculations presented in sections flll Al and lHI Bl 

are a demonstration for the case of numerical relativity 
simulations of coalescing binary compact objects; gravi- 
tational waveform results from different numerical rela- 
tivity codes using different methods and/or different im- 
plementation techniques must be validated in a similar 
way. 

It is instructive to compare the accuracy requirements 
found here with other multiple-orbit binary neutron star 
simulation results, such as those in 0. However, the 
detailed accuracy studies of numerical relativity binary 
simulations in sections IIII Al and IIII Bl are made possi- 
ble by repeating the same simulation (more specifically, 
using the same initial data) many times using a wide 
variety of discretization parameters and outer boundary 
placements, as presented in 0. While the simulations 
presented in [6j used several discretization parameters 
and boundary placements, they were performed for dif- 
fering initial data corresponding to an array of initial 
binary separations. Thus, the detailed studies in sec- 
tions IIII Al and IIII Bl cannot be repeated using the re- 
sults from the multiple orbit neutron star simulations 
presented in [6J. However, a casual inspection of the sim- 
ulations from jfj confirms that the errors induced by the 
boundary are similar to those of 4] . In Fig. [SJ a Simula- 
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tion of initially corotating binary neutron stars from [6( 
is displayed; the coordinate separation of the neutron 
stars is plotted as a function of time over two orbital 
periods. The equations of state used in the simulations 
of Q| and are identical, and the results from plot- 
ted in Fig. |H1 use neutron stars that are 7% more mas- 
sive than those used in .4]- Note that the simulation 
from |(j (labeled "Stable" in Fig.l of reference [6j) dis- 
played in the top panel in Fig. [S] has an initial separation 
of Ti — 9.85 M (separation values in Q are normalized 
by the total baryonic mass M ; we follow this convention 
in Fig. [S] and during the discussion here). After two or- 
bits, the binary separation has increased over 10%, while 
a post-Newtonian point particle simulation using iden- 
tical mass, initial separation, circular orbit initial condi- 
tions, and accurate to order (v/c) 9 in the post-Newtonian 
expansion, predicts that the separation should instead 
decrease by 15% during the first two orbits (the post- 
Newtonian simulation is plotted along side the "Stable" 
simulation from 6] in the top panel of Fig. [SJ. At the 
very least, it is clear that the separation of the neutron 
stars cannot increase, due to the fact that i) the dissi- 
pative effects of gravitational radiation will cause a de- 
crease in the binary separation and ii) while the circu- 
lar orbit initial condition induces a slight eccentricity to 
the orbit of the binary, this initial condition corresponds 
to an apastron (maximum separation) point in the dy- 
namical evolution j 1 31 ] . To compare with the simulations 
analyzed in this paper, the simulations NS-A and NS- 
B from jj] shown in Fig. 01 are reproduced in the lower 
panel of Fig. |SJ Simulation NS-B has a similar outer 
boundary placement as that used in the "Stable" simu- 
lation from 6] , but NS-A has an outer boundary that is 
twice the distance from the binary as compared to NS-B. 
The similarities between the top and bottom panels of 
Fig. IS] suggests that the cause of the unphysical increase 
in the binary separation during the first two orbits of 
the "Stable" simulation from jy] is the close proximity 
of the boundary of the computational domain, and that 
the errors induced by the boundary of the computational 
domain in [4j and |ij are of a similar magnitude. 



IV. CONCLUSIONS 

Using a criterion for gravitational waveform template 
accuracy motivated by matched filtering and parameter 
estimation requirements of modern interferometric grav- 
itational wave detectors, I have calculated the accuracy 
required of numerical relativity simulations of coalescing 
compact binary systems. I have calculated the numerical 
errors of state-of-the-art numerical relativity simulations 
of orbiting binary neutron stars |4| , and I find these errors 
to be several orders of magnitude larger than the allowed 
errors determined from gravitational waveform accuracy 
considerations. Using a post-Newtonian model for the 
truncation errors and boundary errors in the numerical 
simulations of 4j, the computational resources required 
in order that these simulations attain an accuracy needed 
for reliable gravitational wave extraction have been calcu- 
lated. I find that while mesh refinement technology will 
provide an improvement over the unigrid second-order 
accurate simulations of , higher order methods will also 
be required for a robustly accurate numerical relativity 
calculation of multiple-orbit binary coalescence calcula- 
tions on Terabyte-scaled (or smaller) digital computers. 
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